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ABSTRACT 

This paper presents a stochastic approach to the clustering evolution of dark 
matter haloes in the Universe. Haloes, identified by a Press-Schechter-type al- 
gorithm in Lagrangian space, are described in terms of 'counting fields', acting 
as non-linear operators on the underlying Gaussian density fluctuations. By 
ensemble averaging these counting fields, the standard Press-Schechter mass 
function as well as analytic expressions for the halo correlation function and 
corresponding bias factors of linear theory are obtained, extending the recent 
results by Mo and White. The non-linear evolution of our halo population is 
then followed by solving the continuity equation, under the sole hypothesis 
that haloes move by the action of gravity. This leads to an exact and general 
formula for the bias field of dark matter haloes, defined as the local ratio be- 
tween their number density contrast and the mass density fluctuation. Besides 
being a function of position and 'observation' redshift, this random field de- 
pends upon the mass and formation epoch of the objects and is both non- linear 
and non-local. The latter features are expected to leave a detectable imprint 
on the spatial clustering of galaxies, as described, for instance, by statistics like 
the bispectrum and the skewness. Our algorithm may have several interesting 
applications, among which the possibility of generating mock halo catalogues 
from low-resolution N-body simulations. 

Key virords: galaxies: clustering - cosmology: theory - large-scale structure 
of Universe ~ galaxies: formation - galaxies: evolution - galaxies: haloes 

1 INTRODUCTION 

The theory proposed by Press and Schechter (1974, hereafter PS) to obtain the relative abundance of matter con- 
densations in the Universe has strongly influenced all later studies on the statistical properties of dark matter haloes 
and led to a large variety of extensions, improvements and applications. Actually, already in the sixties, Doroshkevich 
(1967) had derived the mass distribution function for 'newly generated cosmic objects', completely analogous to the 
PS one; he had also clearly pointed out the existence of what has been later referred to as cloud-m-cloud problem 
(e.g. Bardeen et al. 1986). The 'Press-Schechter model', which is based on the gravitational instability hypothesis, 
is now considered as one of the cornerstones of the hierarchical scenario for structure formation in the Universe. It 
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shows, in fact, how gravitational instabihty makes more and more massive condensations grow by the aggregation 
of smaller units, only provided the initial density fluctuation field contains enough power on small scales. The main 
drawback of the original PS model is indeed the cloud-in-cloud problem, i.e. the fact that their procedure selects 
bound objects of given mass that can have been already included in larger mass condensations of the same catalog. 
The problem was later solved by several authors (Peacock & Heavens 1990; Bond ot al. 1991; Colo 1991) according 
to the so-called 'excursion set' approach, by calculating the distribution of first-passage 'times' across the collapse 
threshold for suitably defined random walks. Lacey and Cole (1993, 1994) implemented these ideas to study the 
merger rates of virialized haloes in hierarchical models of structure formation. 

An important aspect of the PS model is that, being entirely based on linear theory, suitably extrapolated to the 
collapse time of spherical perturbations, it is, by definition, local in Lagrangian space. While this Lagrangian aspect 
of the theory does not have immediate implications for the study of the mass function of dark matter haloes, it is, 
instead, of crucial importance for their spatial clustering properties. This point was recognized by Cole and Kaiser 
(1989) and, more recently, by Mo and White (1996, hereafter MW), who proposed a bias model for halo clustering 
in Eulerian space, by a suitable extension of the original PS algorithm for the mass function. With their formalism 
MW studied the clustering of dark matter haloes with different formation epochs (sec also Mo, Jing & White 1996). 
The comparison of their theoretical predictions with the spatial distribution of haloes obtained by a friends-of -friends 
group finder and by a spherical overdensity criterion in numerical simulations proved extremely successful. 

These very facts imply that there exists a local version of the PS algorithm providing a mapping between points 
of Lagrangian space and the haloes in embryo which will come into existence at the various epochs. For a given 
realization of the initial density field, the PS mapping is such that, at a fixed redshift z, each Lagrangian point q can 
be assigned to a matter clump of some mass M, identified by a suitable Lagrangian filter, which is collapsing at the 
epoch Zf = z. One can therefore exploit the existence of this mapping to assign a stochastic halo process, our halo 
counting field below, to each point q. This will be the starting point of our analysis. 

What the PS ansatz cannot account for is the fact that the fluid elements are moved apart by gravity, so that 
the halo which the PS mapping assigns to the fluid patch with Lagrangian coordinate q is not going to collapse in 
the same position, i.e. at x = q, but, rather in the Eulerian point x(q, z) — q + S(q, z), with S(q, z) the displacement 
vector, corresponding to the Lagrangian one at the epoch z = z/(q, M) of its collapse. This fact, while not affecting 
in any way the PS result for the mean mass function, as the average halo abundance cannot change by scrambling the 
objects, sensibly modifies their spatial clustering properties. Modeling the latter effect is one of the main purposes of 
the present work. In their derivation of the Eulerian halo bias MW took into some account this problem by allowing 
for the local compression, or expansion, of the volumes where the haloes are located, an effect which is of crucial 
importance for the derivation of the correct halo density contrast. Their derivation, however, is formally fiawed by 
the fact that they only deal with mean halo number densities, so that they are forced to define the bias in terms of 
them. For reasons to be shown below, however, this heuristic treatment can be put on sounder statistical grounds, 
by applying a suitable coarse-graining procedure. 

Of course, the PS model has its own limitations. The comparison of its predictions for the mass function with 
the outputs of N-body simulations (e.g. Efstathiou et al. 1988; Gelb & Bcrtschingcr 1994; Lacey & Cole 1994), while 
surprisingly successful in its general trends, given the simplicity of the assumptions, showed a number of problems. 
Gelb and Bertschinger (1994), for instance, found that the simulated haloes are generally less massive than predicted, 
the reason being that merging does not erase substructure in large haloes as fast as required by the PS recipe. 

There have been many attempts to improve the original PS model. If cosmic structures preferentially formed at 
the peaks of the initial density fluctuation field this would affect their mean mass function (Bardeen et al. 1986; Bond 
1988; Colafrancesco, Lucchin & Matarrese 1989; Peacock & Heavens 1990; Manrique & Salvador-Sole 1995, 1996). 
Bond and Myers (1996) developed a peak-patch picture of cosmic structure formation, according to which virialized 
objects are identified with suitable peaks of the Lagrangian density field. The peak-patch collapse dynamics is then 
followed in terms of the homogeneous ellipsoid model, which allows for the influence of the external tidal field, while 
the Zel'dovich approximation (Zel'dovich 1970) is used for the external peak-patch dynamics. The effects of non- 
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spherical collapse on the shape of the mass distribution were studied by Monaco (1995). Lee and Shandarin (1997) 
analytically derived the mass function of gravitationally bound objects in the frame of the Zel'dovich approximation. 

We prefer here to follow the simple lines of the PS theory to set up the 'initial conditions' for our stochastic 
approach to the evolution of halo clustering. Nevertheless, one should keep in mind that our approach is flexible 
enough to accept many levels of improvement in the treatment of the Lagrangian initial conditions. 

A relevant part of the following analysis will be devoted to the study of the evolution of halo clustering away 
from the linear regime. It turns out that the problem can be solved exactly in terms of the evolved mass density. An 
important result of this analysis is that the general Eulerian bias factor, deflned as the local ratio between the halo 
density contrast and the mass fluctuation field, is both non-linear and non-local. The latter property follows directly 
from our selection criterion of candidate haloes out of the linear density field. 

Our algorithm can also be seen as a specific example of a bias model which is local in Lagrangian space. This is 
expected to have relevant consequences on galaxy clustering. Because of this local Lagrangian character, our model 
strongly differs from the local Eulerian bias prescription applied by Pry and Gaztafiaga (1993) to the analysis of 
the hierarchical correlation functions. A simple test of our theory can be obtained by analyzing the behaviour of the 
bispectrum (or the skcwness), whose shape (scale) dependence will be shown to be directly sensitive to the assumption 
of local bias in Lagrangian vs. Eulerian space. 

Our results for the evolved halo distribution generally allow to study their statistical properties at the required 
level of non-linearity, and could be further implemented to generate mock halo catalogues starting from low-resolution 
numerical simulations of the dissipationless matter component. These results have important implications for the 
study of the redshift evolution of galaxy clustering, a problem made of compelling relevance by the growing body of 
observational data at high-redshift which are being produced by the new generation of large telescopes. A general 
study of this problem has been recently performed by Peacock (1997) and Matarrese et al. (1997); the latter pointed 
out that knowledge of the evolution of the effective bias for the various classes of objects is a key ingredient in the 
comparison of theoretical scenarios of structure formation with observational data on clustering at high redshift. 
Kauffmarm, Nusser and Stcinmctz (1997) used both semi-analytical methods and N-body techniques to study the 
physical origin of bias in galaxies of different luminosity and morphology. 

The plan of the paper is as follows. In Section 2 we define our halo counting field, within the linear approximation, 
both in the Lagrangian and Eulerian context. The non-linear evolution of the halo clustering is analyzed in Section 
3, where we also compute the bispectrum and skewness of the evolved halo distribution. Section 4 contains a general 
discussion of our results and some conclusions. 



2 STOCHASTIC APPROACH TO HALO COUNTING AND CLUSTERING 

2.1 Basic tools and notation 

Let us assume that the mass density contrast e(q), linearly extrapolated to the present time, is a statistically homo- 
geneous and isotropic Gaussian random field completely determined by its power-spectrum P{k). Here q represents a 
comoving Lagrangian coordinate. A smoothed version of the field e(q) is obtained by convolving it with a rotationally 
invariant filter Wii(g), containing a resolution scale R, with associated mass M p^B?, pi, being the background 
mean density at = 0, 

e«(q) = j rfq' Wr{\^ - q'l) e(q') = €^(q) . (1) 

The smoothed field is also Gaussian with one-point distribution function Ga^i^M) ~ (2"" )~^''^ ^^p(~^m/2''"m)' 
where cr^ denotes the variance of e^, = (e^) = (27r^)~^ dkk^ P{k)W{kR)^ . The symbol W{kR) indicates 
the Fourier transform of the filter function. In the following, we will often be concerned with the joint distribution of 
the fields eMj(q) and e„^(q). The two-point correlation function of the linear density contrast smoothed on the scale 
Ri and R2 is 
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1 . . 

Ci2(g) = {e„^(qi)6„^(q2)) = 2^y^ dkk^P{k)W{kR,)W(kR2)Mkq) , (2) 

where g = — qjl and jo (a;) is the spherical Bessel function of order zero. We term (T12 the value assumed by ^12 
in the limit g — » 0. 

The properties of the filtered quantities clearly depend upon the choice of the window function. For instance, the 
relation between the mass enclosed by a top-hat filter Wfl(g) — 30{R — q)/4-KR^ (where Q{x) is the Heaviside step 
function) is the standard M(R) — A-KpbR? /?>■ Instead, for a Gaussian window, WR{q) — (2nR^)~^^^ exp{—q^ /2R^), 
the enclosed mass is M{R) = {2Tvf^^ptR^. These two masses coincide for Re ~ 0.64 i?TH (Bardeen et al. 1986). 

In the literature, the sharp top-hat filtering has been alternatively adopted in Fourier space, Wii{k) = 0(fc_R — k), 
where ka = 1/R and k — |k|. The most remarkable property of this filter is that each decrease of the smoothing 
radius adds up a new set of uncorrelated modes (Bardeen et al. 1986; Bond et al. 1991; Lacey & Cole 1993). This also 
implies that, for example, the correlation function in eq.(^ simplifies to ^12 = ^ii, whenever fenj < kn^ ; consequently, 
o"i2 = (Til = o"i. In practice, the information is always erased below the largest of the two smoothing lengths. This 
property will be particularly useful in the next sections. For this 'sharp k-space' filter, the main difficulty is how to 
associate a mass M{R) to the cutoff wavenumber kR. Lacey and Cole (1993) give the expression M{R) = Qtv^ p^k^"^ , 
which coincides with the mass within a top-hat filter if one takes kji = 2.42/ Rth ■ 

In the next section we introduce the halo counting random fields that allow a fully stochastic description of the 
biased haloes distribution. To illustrate how our formalism works, we first show how to derive the PS mass function 
by performing a simple averaging of our stochastic counts. 



2.2 Lagrangian mass function: Press-Schechter theory 

Press and Schechter proposed a simple model to compute the comoving number density of collapsed haloes directly 
from the statistical properties of the linear density field, assumed to be Gaussian. According to the PS theory, a patch 
of fluid is part of a collapsed region of scale larger than M{R) if the value of the smoothed linear density contrast on 
that scale exceeds a suitable threshold t / . The idea is to use a global threshold in order to mimic non-linear dynamical 
effects ending up with halo collapse and virialization. An exact value for tf can be obtained by describing the evolution 
of the density perturbations according to the spherical top-hat model. In this case, a fluctuation of amplitude e will 
collapse at the redshift Zf such that e(q) = tf = 5c/D{zf), where D{z) denotes the linear growth factor of density 
perturbations normalized as D{0) = 1. In the Einstein-de Sitter universe and during the matter dominated era the 
critical value Sc does not depend on any cosmological parameter and is given by 5c = 3(1271)^/^20 ~ 1.686, while, 
for general non-flat geometries, its value shows a weak dependence on the density parameter Q, the cosmological 
constant A and the Hubble constant H (e.g. Lacey & Cole 1993), thus on redshift. In a flat universe with vanishing 
cosmological constant D{z) = {1 + z)~^; explicit expressions for the linear growth factor are given in Appendix A for 
general Friedmann models. 

A local version of the PS approach can be built up in terms of stochastic counting operators acting on the 
underlying Gaussian density held, as follows. The number of haloes per unit mass, contained in the uiut comoving 
volume centered in q, identifled by the collapse-threshold tf{zf), is described as a density held of a point process by 

AA,f (qlM, tf) = -2 II [e„(q) - tf] . (3) 

Note that the quantity J\f^{q\M, tf) is non-zero only when the flltered density contrast in q upcrosses, or downcrosses, 
the threshold tf,hy varying the smoothing length R (or the corresponding mass M). The factor of 2, appearing in the 
expression of A/'^(q|M, t/), is needed in order to obtain the right normalization of the mass function, in which case 
it has been shown to be intimately related to the solution of the cloud-in-cloud problem (Peacock & Heavens 1990; 
Bond et al. 1991; Cole 1991), at least for sharp k-space flltering. At this level, our description should be thought of as 
a sort of differential version of Kaiser's bias model (Kaiser 1984), that defines a population of objects with the right 
average halo abundance and their related clustering properties, rather than a detailed modeling of how structures 
form from the primordial density field. In a forthcoming paper, however, we will show that the present approach is 
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fully consistent with a rigorous treatment of the cloud-in-cloud problem (Porciani et al. 1997) . In that approach halo 
correlations will be obtained from pairs of first-upcrossing 'times' for spatially correlated random walks above the 
collapse threshold tf. 

It can be seen from equation (^) that a population of haloes is uniquely specified by the two parameters M 
and tf. In the standard PS formulation tf is interpreted as a sort of time variable, related to the formation redshift 
Zf, which decreases with real time, as every halo continuously accretes matter. In this sense one can say that, for a 
continuous density field with infinite mass resolution, each halo disappears as soon as it forms to originate another 
halo of larger mass. 

Alternatively, instead of considering tf as a time variable, one can use it simply as a label attached to each halo. 
The haloes so labelled can be thought as keeping their identity during the subsequent evolution at any observation 
redshift z. This is not in contrast with the fact that in the real Universe dark matter haloes undergo merging at 
some finite rate (e.g. Lacey & Cole, 1993, 1994). Within such a picture, in fact, the physical processes of accretion 
and merging reduce to the trivial statement that haloes identified by a given threshold are necessarily included in 
catalogues of lower threshold, so that, in the limit of infinite mass resolution, only haloes with Zf — z would actually 
survive. Nevertheless, keeping Zf distinct from z may have several advantages, among which the possibility of allowing 
for a more realistic description of galaxy and cluster formation inside haloes, for both the evolution of the luminosity 
function (Cavaliere, Colafrancesco & Menci 1993; Manrique & Salvador-Sole 1996) and of the galaxy bias (e.g. MW; 
Kauffmann et al. 1997). Let us stress, however, that we are not addressing here the issue of galaxy or cluster merging: 
our method is completely general in this respect and allows to span all possible models, from the instantaneous 
merging hypothesis {zf — z) to the case of no merging at all (z/ fixed for changing z < Zf). 

In what follows, therefore, we will assume that we can deal with the halo population specified by the formation 
threshold tf at any redshift z. Only in this sense we will say that we 'ignore' the effects of merging in our description: 
merging can be exactly recovered at any step, and with any assumed mass resolution, as the relation between Zf and 
z. To implement this idea it is enough to scale appropriately the argument of the Heaviside function in eq.(^, which 
can be recast in the form 

J^,^ici,z\M,Zf) = ~2^^e[eJci,z)-Sfiz,Zf)] , (4) 

where ej^j{q,z) = D(2:)e(q) and Sf{z,Zf) = ScD{z)/ D{zf). It can be easily seen that the ensemble average of the 
counting field A//f (q, z\M, Zf) corresponds to the PS mass function 

{Ntici, z\M, Zf)) dM = {z\M, Zf) dM , (5) 
where 



n (z\M Zf) dM = — ^ "^^(^'^-f^ exD 



S']iz,zf) 



dM 



dM . (6) 



Note that we emphasized the z-dependence of the comoving mass function, though it is straightforward to verify that 
the value of npg{z\M,Zf) does not change with z. In fact, since we are ignoring the effects of merging, once a class 
of haloes has been identified, their mean comoving density keeps constant in time. Thus, as far as the mass function 
is concerned, the introduction of the observation redshift z is somewhat more formal than physical. However, this 
distinction will be far more significant in the next sections, where, in order to compute the halo-to-mass bias factor, 
we will relate the Lagrangian distribution of a population of haloes selected at Zf to the mass density fluctuation 
field linearly extrapolated to the redshift z. Models of galaxy formation which assume that galaxies form at a given 
redshift Zf with some initial bias factor and that their subsequent motion is purely caused by gravity (e. g. Dekel 
1986; Dekel & Rees 1987; Nusser & Davis 1994; Fry 1996) can be easily accommodated into this scheme. 
To conclude this section, let us consider the integral stochastic process 

dM'M'Ht{ci,z\M',Zf) ^2p,Q[e,,{q,,z)-5f{z,Zf)] , (7) 
representing the fraction of mass, in the unit Lagrangian comoving volume centered in q, which at redshift Zf has 
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formed haloes more massive than M. This coincides with the original Kaiser bias model (Kaiser 1984) up to the 
multiplicative factor 2pf,, which is irrelevant for calculating correlation functions. 



2.3 Conditional Lagrangian mass function 



The PS theory reviewed in the previous section describes the overall distribution of halo masses in a homogeneous 
universe of mean density pb- However, of cosmological interest is also, for instance, the estimate of the halo distribution 
within rich or poor environments (which can be related to the galaxy number enhancement per unit mass in rich 
clusters or in voids), thus justifying the investigation of the distribution of halo masses conditioned to lie within a 
larger uncoUapsed container of given density. The conditional mass function has been studied by several authors (e.g. 
Bond et al. 1991; Bower 1991; Lacey & Cole 1993). 

We extend here the approach introduced in the previous section in order to derive the conditional mass function. 
Specifically, we calculate the comoving mass function, in the mass range M, M + dM, for objects contained in a large 
region of dimension Ro, corresponding to a mass Mo, with local density contrast to = ^j^j^ ■ We will require eo <^ Sf 
and Ro 2> R, to ensure that the container is not collapsed yet by the epoch zj and that it encloses a non- negligible 
population of objects. 

In order to mimic these environmental effects, we modify the halo counting field according to 



Afh{ci,z\Ad,Zf\Mo,eo) ^ -^^^Mi^l^z) - Sf(z,Zf)] 5„[ej,,j^^(q, z) - to] , 



(8) 



where ^^(q) denotes the Dirac delta function and A'^o = (^j-, [e„^(q, 2) — eo]) is the normalization constant. Here the 



scalar eo indicates the value of the random field e^^-^ (q, z). Taking the ensemble average (and using the cross- variance 
aij for a sharp k-space filter) one eventually obtains 



(AAfe (q,2|Af, Zf\Mo,eo)} dM = nps{z\M, Zf\Mo,eo)dM 



(9) 



where 

np^{z\M, Zf\Mo,eo 



dM : 



Pb Sf{z,Zf) ■ 



M [ol(z) - ol{z)f/^ \ 2 [aliz) - al{z)] 

This straightforward calculation shows how to obtain results already known in the literature by simply starting from 
the random field in eq.(^: averaging that halo counting field leads to the expected conditional mass function. But 
not only: unlike previous treatments, once the halo counting field has been consistently defined, other statistics, like 
the two-point halo correlation function, can be calculated. We will carry out this program in the next section. 



dM 



dM 



(10) 



2.4 Lagrangian clustering: halo-to-mass bias from correlations 

In this section we will compute the halo-halo correlation function which coincides with the correlation function of our 
random counting field. Specifically, we will calculate the Lagrangian halo correlation function from the Lagrangian 
counting field A/jf (q, z\M, Zf). By definition, the correlation function of this stochastic process is given by 

{Khl,z\Ml,5f{z,Zl)]M^[q2,z\M2,5f{z,Z2)]) 



^hh{q,z\Mi,zi;M2,Z2) = 



1 , 



{J^^[q„z\Ml,5f{z,Z,)]){J^,^[ci„z\M2,Sf{z,Z2)]) 

where q = |qj^ — qj]. Performing the ensemble average over the Gaussian fields e^^^ (q) and e^j^^ (q), we obtain 



(11) 



{M^[ci„z\M,,5f{z,Z,)]M,f[ci„z\M2,Sf{z,Z2)]) = 



d d 



M1M2 dMi dM2 
where G2 (01,02) denotes the bivariate Gaussian distribution 



dofi da2 G2(qi, Q2) , (12) 



G2(Q1, Q2) = 



exp 



— 2u) 

j| (Tl (T2 



(13) 



with normalized correlation ioici) = Ci2(g)/o"Mj ''"m^ and Oi = D{z)o^j, . 
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The full exact expression for the halo-halo correlation function can be obtained after an incredibly long algebraic 
computation. We report here only the final expression. Defining Sfi = 5f{z, Zi), we have 



l + ^L{<l,z\Mi,zi;M2,Z2) 



da2 



+ 



+ 



{do 

\ dMl dM2 ' 5f2{l-Lu'^) V cTi 
'5/2 5/1 \ du! da2 



V CTi (72 / 



5/2 \ dai dio 



5/1 (1-^2) 

(7l(72 
5/15/2(1-^^2)2 



\ (72 (7i 



+ 



CTl(72 



dMl dM2 



dMl dM2 5/i5/2 dMidM2 
5?i 



a;(l-a;^) + (l+^^)^^ - 

(7l (72 



+ 



9cj 9cj 



aMi 9M2 



exp 



5^ (ga 

erf (7| 



(71 (72 



2(l-a;2 



ci(7i C?(72 



/ dai 
[dMl 



dM2 



(14) 



This expression can be easily shown to be independent of the observation redshift z. A remark is now appropriate. 
Our formalism describes the halo distribution as a discrete point process. However, actual haloes are extended in 
size. This is why, as also seen in numerical simulations, for separation smaller than the typical Lagrangian radius of 
the halo, the correlation function abruptly reaches the value —1: a sort of 'exclusion principle' for extended haloes. 
Thus, we expect that the correlation function in eq.(|lj) can be a reliable description of halo clustering only for 
q Max(f?i, 7?2). Another point concerns the use of finite mass resolution as in A'^-body simulations. The proper 
analytical correlation to compare with in that case is the integral of Chh ^ps (-^i) (-^2) over the specified mass 
interval, appropriately normalized. 

Since the action of the window functions on the correlations is negligible for lags q much larger than the smoothing 
lengths, q ^ Ri and g S> -R2, for the normalized correlation we obtain a;((/) ~ Cm (g) /(7jv,^ cr^j^ [where ^miq) is the 
linear mass autocorrelation function] and, eventually, for the halo correlation 



Lh{q,z\Mi,zi;M2,Z2 

Explicitly, the first two bias parameters read 
Sf{z,zf) 1 _D(zf) 



h'i{z\Mi,zi)h'({z\M2,Z2)^„,{q,z) + -bf^2{AMi,Zi)hii{z\M2,Z2)iUq,z) + - 



K{z\M,Zf) = 



h^2(AM,Zf) = 



al{z) 5f(z,zf) D{z) 



5j{z,Zf) 



(^) 



51 



D{z) 



D{zf) 



(15) 



(16) 



(17) 



These expressions for the bias factors generalize, in a sense, those concerning the clustering properties of dark matter 
haloes in Lagrangian space obtained by MW and Mo et al. (1996), with the relevant difference that we have obtained 
the bias factor from the behaviour of the halo two-point correlation function. Also relevant is the fact that, unlike 
MW, we obtained our Lagrangian correlation function without introducing any background scale Ro, which allows 
to extend its validity down to spatial separation comparable to the halo size R <ti Ro. A calculation of the leading 
behaviour of the correlation deriving from equations ( |ll| ) and ( |l^ has been already attempted by Kashlinsky (1987) 
who, however, missed the contributions originated by the differentiation of cu with respect to AIi and M2, thereby 
obtaining an incomplete expression for bi . 

The halo correlation function in Lagrangian space, S^f^f^ from eq.(p^, with Mi — M2 — M and 2:1 = 22 = 2/, 
is shown in Figure 1 for two scale-free power-spectra, P{k) oc fc", with spectral index n — —2 and n — —1, in the 
Einstein-de Sitter case. The two-point function is calculated for various halo masses in units of the characteristic mass, 
M, , defined so that (7m* =tj — 5c/D{zf), with top-hat filtering; the spatial dependence is shown as a function of 
the scaling variable q/R, which eliminates any residual redshift dependence. Also shown is the mass autocorrelation 
function and an estimate of the Lagrangian halo two-point function obtained as {bi)^$,m, for M 7^ M,, and 
(^2)^Cm/2, for M = M«, as in this case hi = 0. Note that such an estimate of always provides an accurate fit to 



* We are adopting here the MW definition of Af, , which, although differing from the standard PS one, cri\/* = tf/\/2, is more 
convenient for our present purposes. 
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1 1 — I — I I I I 1 1 — r3 



M/M,= l 



~i 1 — I — I I I I 1 1 — r3 




M/M.= l/8 



_l I I I I L_ 



n — I — I I I I I 1 1 — r 




Figure 1. The exact Lagrangian halo correlation function in an Einstein-de Sitter universe {solid lines) is shown for scale-free 
models with spectral index n = —1 and n — —2, and for various masses. In each panel we set M\ = M2 = M and z\ = Z2- 
Results are plotted in terms of the scaling variables M/M-t and q/R (with R the top-hat radius corresponding to the halo 
mass Af), which makes the resulting curves redsliift independent. For comparison, the linear mass autocorrelation function 
smoothed on the halo scale is also shown with long-dashed lines. The short-dashed lines represent the linear bias approximation 
for the halo correlations: (fe^)^^m. In the central panels, where = 0, the dot-dashed lines show, instead, the second-order 
approximation for . Each column contains panels that refer to the same mass variance (c^/ij = 1/4, 1, 4) and so to the 
same Lagrangian bias factors. Notice that, for separation a few times the halo size the first non- vanishing term of eq.(15) always 
gives an accurate approximation to the exact halo correlation. This implies that, for M* objects, « (''2 )^€m/2- 



its analytical expression for separation a few times larger than the halo size. The characteristic behaviour of the halo 

correlation function for M = Af, , where the linear bias vanishes, is actually a peculiarity of the Lagrangian case (see 
also M W) . As we will see below, the Eulerian halo correlation function does not show such a drastic change of slope 
in the same mass range. 



2.5 Peak-background split 

In the previous paragraphs we computed number densities and correlation functions of haloes in Lagrangian space. 
However, after their identification, these haloes in embryo move following the gravitational field, modifying their 
original spatial distribution. One issue to address is how, for instance, the conditional halo number density per unit 
mass changes £is a consequence of gravitational evolution. Furthermore, of interest is to quantify the evolution of 
clustering in terms of the halo correlation functions, or in terms of the halo-to-mass bias. Both problems can be dealt 
with by defining Eulerian halo counting fields, in the same spirit as we did for the Lagrangian case. 

Essentially, our approach to the clustering evolution can be based on a generalization of the so-called peak- 
background split, first proposed by Bardeen et al. (1986), which basically consists in splitting the mass perturbations 
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in fine-grained (peak) and coarse-grained (background) components^. The underlying idea is to ascribe the coUapse of 
objects on small scales to the high-frequency modes of the density field, while the action of large-scale structures on 
these non-linear condensations is due to the remaining low-frequency modes. At the linear level the resulting effect 
of these long wavelengths is simply modeled as a shift of the local background density. 

In the spirit of the peak-background split, we define the linear density field smoothed on a given scale em as 
consisting of two complementary and superimposed components, namely ej(/ = etg-l-epfe. Adopting as window function 
the sharp fc-space filter, we define as 'background' component the density contrast smoothed on the scale Ro = 1/fco 

e,,(q, z)^ j ^ ?(k, z) e(fco - k) e^"^ . (18) 

The 'peak' component is instead obtained by smoothing the mass density fiuctuation with the band-pass filter 
©(fcjvf — k) 0{k — ko), namely 



epfe(q, z) = J ?(k, z) e{kM - k) e(fe - feo) e"^ "^ , (19) 

where kM ~ 1/Rm, with M oc pbR%i and Ma oc pbRo the masses enclosed by the two filters. So, the peak component 
contains only modes with wavenumber in the interval [ko, kAi]- Note that in the linear regime, with Gaussian initial 
conditions, the peak and background components are statistically independent, i.e., 

{epk{cij^,z) ebg {ci2, z)) =0 , (20) 

for, by construction, the two fields do not share any common Fourier mode. To summarize: provided the collapsed 
object is described according to the spherical model, as in the PS theory, the peak field e^^ (q, z) can be thought as 
evolving in a local environment with effective mean density pb[l + e(,g(q, 2)]. This implies that the collapse condition 
can be written as epfe(q, z) = 5f{z, zf) — ebg{q_, z). 



2.6 Eulerian halo counting field and bias 

The previous analysis shows how the PS and the conditional Lagrangian mass functions can be obtained by averaging 
properly defined halo counting random fields. It is thus legitimate to explore the possibility of building up analogous 
counting processes in the Eulerian world. Our approach here will be based on the peak-background split technique 
described above. 

Let us define the Eulerian counting field of haloes collapsed at redshift 2/ and observed at 2 as 

J^^iq, z\M, Zf) = [1 + e,,(q, 2)] AA^^(q, z\M, Zf) = [l + e5,(q, 2)] ^6 [^^^(q, 2) - {Sf{z, Zf) - e,,(q, 2))] .(21) 

The watchful reader might wonder about our use of the Lagrangian variable q within the Eulerian framework, however, 
in linear theory, x = q. Once again, the redshift z must be thought of as the redshift the sampled objects have at 
the epoch of observation. Noteworthy, eq.(^) is fully consistent with the analysis in Cole and Kaiser (1989). Most 
importantly, our treatment allows for a local description. Let us stress here that the factor (1 + etg), connecting the 
Eulerian to the Lagrangian counting field, simply comes from mass conservation in Eulerian space [see also Section 
3.1 and, in particular, eq.(^)]; this point has been discussed in more detail by Kofman et al. (1994). 
Now consider the integral stochastic process 

poo 

/ dM'M'U^{ci,z\M',Zf)^2pb[l + ebgici,z)]e[eM{<i,z)-5f{z,Zf)] ; (22) 

J M 

this represents the fraction of mass, in the unit Eulerian comoving volume centered in q, which at redshift 2/ will 
form haloes more massive than M. For Mo —> M, ebg and the above relation coincides (up to the usual fudge 

factor of 2, having no effect on correlations) with the weighted bias model of Catelan et al. (1994). An extended 
version of this scheme, called 'censoring bias', has been recently proposed by Mann, Peacock and Heavens (1997). 
Thus, the weighted bias is just the Eulerian version, within linear theory, of Kaiser (1984) bias model. 

t We are here making a rather liberal use of the word 'peak', to mean the fine-grained component of the linear density field. 
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Of course, further specifications could be added to our Eulerian counting field. For instance, we might ask that 
the background scale has not yet collapsed by the epoch zf, in such a case we should multiply the above stochastic 
process by the factor 0[(5/ (z, z/) — ei,g(q, z)]. Extra details of this kind would however make negligible changes to our 
final results, provided a'^^ ^ ctq. 

Like in the Lagrangian case, to calculate the mean halo number density per unit mass, one needs to ensemble 
average A/^f (q, z|M, z/). Let us analyze this operation in more detail. Because of the way the Eulerian counting 
process has been defined, it is clear that Af^^ {q, z\M, Zf) depends on two random fields, specifically etg and tp^. So, 
the ensemble average {JV^) can be interpreted as a double average over these fields, i.e. (Aff^) = { {■N'h)tpk)ibg- The 
statistics of the field A/'/f can be described in terms of n-th order correlation functions, { (A/^f (q^) • • ■ Afh iSAn)) tpk) <ibg ■ 
The exact calculation of these quantities is rather difficult. However, because of the short-scale coherence of the peak 
field, implied by the 'infrared' cutofi' at ko, its covariance (epfe(qi) fp* (q; +'"))epfc vanishes whenever r ^ Ro, so that we 
can simplify the general halo correlations above as { (A^f (qj ■ • ■ A/'f (q^))^^ J^.^ « { (AA/f (qi))^^^ ■ ■ ■ {A/'^(q„))Ep J^.^ , 
provided we consider sets of points q,, i = 1, . . . , n, with relative separation r^j = Iq^ — q^l ^ Ro- Therefore, with 
the purpose of calculating the mean Eulerian halo number density per unit mass and Eulerian halo correlations, we 
can make the replacement My^ — > {■N'h)t-pk = ^ with only negligible loss of accuracy. According to the definition 



of A/|j in eq.(|2l[), the latter ensemble average gives 



N^{ci,z\M,Zf) 



Pb 



2-K M 



[l + e6<,(q,2)] 



5f{z,Zf) - eftg(q,2) 



exp 



\5fi.z,Zf) - ebg(q,z)]^ 
2\ol{z) ~ oi(z)] 



dM 



(23) 



which, having averaged over the fine-grained mass fluctuations, represents a sort of coarse-grained halo counting 
field. Notice that the fine-grained ensemble average has replaced the original step-function operator of eq.(|2l|) by a 
smoother function, which can then be consistently expanded in series of the background field, as shown below. 

Let us stress that the expression in eq.(|2^) is just the Eulerian analog of eq.(lO) in MW, but the field ehg is here 
a true random field, and so is the process A'^^. The knowledge of allows us to define the Eulerian halo number 
density fluctuation as 

iVf (q, z\M, zj) - (iVf (q, z\M, Zf)), 



Sh i<i,z\M, Zf) 



{N^{<l,z\M,Zf)),, 



h (q,z|M,z/)e6g(q,z) 



(24) 



where we introduced the Eulerian 'bias field' (q, z\M, zf). The second equality in the above equation does not 
mean that the Eulerian fluctuation field is proportional to the background density field etg. In fact, b^ in general 
depends upon e{,g itself. Its functional dependence can be understood by expanding Nl^{q, z\M, z/) in powers of e{,g 
to obtain 

Sh(q,z\M,Zf) = bf(z\M,Zf)ebg{q,z) + ^bf(z\M,Zf)elg{q,z) + ... 

= [l+bi{z\M,Zf)]etg{q,z) + ^[bi;iz\M,Zf) + 2b^{z\M,Zf)]el{q,z) + ... , (25) 

where, for cr^^ ^ CTq, the flrst and second-order Lagrangian bias factors fcf and 62 are those of eq.(p^) and eq.(p^, 
respectively. Accounting for the transformation from the Lagrangian to the Eulerian distribution (e. g. Kofmann et 
al. 1992), one has {A^(q, z|M, z/))^^^ = np^{z\A4, Zf). It can be useful to give explicit expressions for the first two 
Eulerian bias parameters of linear theory 



6f(z|M,z,) = l + ^ 



D{z) 



biiz\M,Zf) 



D{zYal 



- 3 



+ 



D{z) 



D{zjYoI 



(26) 



(27) 



The set of linear theory Eulerian bias factors bf{z) can be obtained from the Lagrangian ones according to the 
general rule 



■ebti + bf 



(28) 



The same method can be applied to the Lagrangian expression, in the sense that we can obtain, similarly. 
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Nt{ci,z\M,Zf) 



Pb Sf{z,Zf) - ebg{ci,z) 



One has, exactly, {N^^ {(\, z\M, Zf))^^^ = {JV^ {q, z\M, zj-)) = rip^ (z| Af, z/). By expanding the coarse-grained La- 
grangian counting field N^{q, z\M, Zf) we can define Lagrangian bias factors at any order. For a^^ ^ these turn 
out to be identical to those deriving from the expansion of the halo correlation in Lagrangian space, eq.(|l^). This 
suggests, however, that these bias factors can be used to describe halo clustering on distances r > R, without any 
further restriction introduced by the background scale Ro. 

The very fact that, for practical purposes, we can replace the exact operator JV^ by the locally averaged one 
demonstrates that the MW treatment can be made self-consistent, provided their small-scale density field is replaced 
by the peak field, and the value of the threshold is modified accordingly. Most importantly, our local averaging 
procedure implies that, up to the scale Ro, we are indeed correctly accounting for the cloud-in-cloud problem. This is 
because at each point q, characterized by a random value of the background field e6g(q), the coarse-grained stochastic 
process {q, z\M, zj) (and its Lagrangian equivalent) actually represents the local mean mass function, for which 
the cloud-in-cloud problem is exactly solved in terms of first passage 'times' across the local barrier 5f{z, z/) — e6g(q, z), 
with initial condition epfe(q, z) = at _R = i?o. Therefore, to the aim of calculating correlations on lags r ^ Ro, we 
can safely state that our coarse-grained halo counting fields are unaffected by the cloud-in-cloud problem. 

The shift by 1 of the linear bias factor, here implied by the transformation from the Lagrangian to the Eulerian 
world, was also noticed in the weighted bias approach by Catelan et al. [1994; their eq.(21)], where an underlying 
Lognormal distribution was assumed to avoid negative-mass events. 

The above expression for bf {z\M, zf) coincides with the formula by MW [their eq.(20)], who, however, only 
presented it for z = 0. As noticed by MW, an important feature of this linear bias is that it predicts that large-mass 
objects (actually those characterized by < tf) are biased with respect to the mass [bf > 1), while small-mass ones 
(ctj,j > tf) are anti-biased (bf < 1). Haloes with mass close to the characteristic one, M*, have non-vanishing linear 
bias, unlike the Lagrangian case. As we will see in Section 3.1, this one-to-one classification of biased and anti-biased 
objects according to their mass is no longer valid in the non-linear regime, as the shear field at the Lagrangian location 
of the halo also contributes to the determination of its Eulerian bias factor. 

The effect of merging can be easily accommodated into this scheme. In the real Universe, haloes undergo merging 
at some finite rate, which can be suitably modeled (e.g. Lacey & Cole 1993). As mentioned above, in the simple PS 
theory such a rate is actually infinite, for infinite mass resolution, implying that only haloes 'just formed' can survive, 
so that Zf — z. So, if one gives up singling out the individuality of haloes selected at different threshold, i.e. with 
different formation redshifts Zf > z, one immediately obtains (e.g. Matarrese et al. 1997) 



dM 



(29) 



bf{z\M) = 1 + 



Dizr 



which implies a quadratic redshift dependence in the Einstein-de Sitter universe. 



bf{z\M) = 1 + 



5c{l 



1 



(30) 



(31) 



The latter form coincides with the result by Cole and Kaiser (1989) [their eq.(6)], who however define the bias factor 
of haloes at redshift z with respect to the mass fluctuation at the present time, which then scales the latter expression 
by a factor (1 -I- z)~^. 

On the other hand, for fixed Zf and varying z, i.e. for objects which survived till the epoch z after their birth at 
Zf, the Eulerian bias of eq.(^) gets a completely different evolution, namely 

Dizf) 



bf(z\M) = 1 + 



D{z) 



b?{zf\M) 



which implies a linear redshift dependence in the Einstein-de Sitter case, 
1 



bf{z\M) = l + 



6f(2/|A/)-l 



(32) 



(33) 



The latter form coincides with that obtained by Dekel (1986), Dekel and Rees (1987), Nusser and Davis (1994) and 
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Fry (1996). This relation can be relevant for galaxies which were conserved in number after their formation, i.e. that 
maintained their individuality even after their hosting haloes merged. 

It is trivial, at this point, to obtain the Eulerian halo-halo correlation function within our approximations. For 
lags r ^ Ro, one has 

Cf';,(r,^|Mi,2i;M2,Z2) =bf(2!Mi,zi)6f(^lA/2,^2)Cm(r,2) + i6f(z|Mi,zi)6f(^|M2,Z2)d(r,z) + --- . (34) 

The main limitation of this formula, however, is that it only provides a link between the Eulerian halo correlation 
function and that of the mass within linear theory. What one would really need, instead, is a similar relation in the 
fully non-linear regime. This problem will be solved in the next section. 



3 HALO COUNTING AND NON-LINEAR DYNAMICS: EULERIAN DESCRIPTION 

One can derive a general expression for the Eulerian halo-to-mass bias by integrating the continuity equations for the 
mass and for the halo number density, assuming that haloes move according to the velocity field determined by the 
matter. The Lagrangian analysis carried out in the previous section is crucial to the present purposes, since it allows 
for the natural initial conditions necessary to integrate the Eulerian equations. As we will show below, the Eulerian 
halo-to-mass bias obtained in such a way holds for any cosmology and in any dynamical regime. This turns out to 
be a remarkable generalization of the biasing proposed by Cole and Kaiser (1989) and MW. 



3.1 Eulerian bias from dynamical fluid equations 

Let us consider the mass density fluctuation field (5(x, t[z)) = (5(x, z) which obeys the mass conservation equation 

^ = -(l+5)V-v, (35) 

where r is the conformal time of the background cosmology and the differential operator d/dr = d/dr + v ■ V is the 
convective derivative. The peculiar velocity field v = djc/dr satisfies the Euler equation dv/dr + [a' /a)w — —\/<j)g, 
where a is the expansion factor and a prime denotes differentiation with respect to r. For later convenience, let us 
also define the scaled peculiar velocity u = dx/dZ) = v/D'. The peculiar gravitational potential (j)g is determined by 
the matter distribution via the cosmological Poisson equation V^4>g — iivGa'^ pb{T)S, where pi,(r) is the background 
mean density at time r. If we assume that our halo population of mass M and formation redshift 2/ is conserved in 
time, and evolves exclusively under the influence of gravity, its number density fluctuation (5;i(x, z) = 5h{:K, z\Al, Zf) 
has to satisfy the continuity equation (e.g. Fry 1996) 

^ = -(l+<5.)V-v, (36) 



from which, eliminating the expansion scalar V ■ v, we obtain 

d\n{l + Sh) _ dln{l + S) 
dr dr 

This equation can be integrated exactly in terms of Lagrangian quantities, and the solution reads 



(37) 



1 + (x, z)=[l + 5h(q)] [1 + <5(x, z)] (38) 

(see also the discussion in Peacock & Dodds 1994), where q is the Lagrangian position corresponding to the Eulerian 
one via x(q, 2) = q + S(q, 2), with S(q, 2) the displacement vector. In eq.(p8[), by (5h(q) — 5ft(qjM, 2/) we mean the 
Lagrangian halo density fluctuation, whereas, for simplicity, we assumed that lim2_»oo <5[x(q, 2), 2] = 5{q) = 0, i.e. 
that the mass was initially uniformly distributed (this amounts to taking purely growing- mode initial perturbations). 
Defining the Eulerian halo bias field through 

(5h(x, 2) = b^(x, 2) (5(x, 2) , (39) 
we end up with the exact relation 
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b^(x,z) = l + i±^i^5,(q). (40) 
d(x, z) 

The key problem now is how to calculate the field (5h(q). We cannot simply take the Lagrangian halo distribution as 
<5h(q) = b^(q)5(q), because (5(q) = 0; thus, we are forced to adopt some limiting procedure. To specify the Lagrangian 
halo distribution, we can take advantage of the results of Section 2. By definition, the Lagrangian distribution of 
nascent haloes of mass M and formation epoch Zf is given by 



Sh{ci\M, Zf) 



lim h'^ {q, z\M , Zf) etgiq, z) = b^ (qjM, 2/) eo(q) 



(41) 



where bo(q|M, «/) = b^(q, 2 = 0|M, z/) is the Lagrangian halo bias field. Once again, let us stress that the second 
equality in the latter equation does not mean at all that 5/i(q) is proportional to eo(q). In fact, b^ is in general a 
functional of the background density field. To understand the above equation, one has to remember that, at sufficiently 
early times, the expression for the Eulerian bias field obtained in linear theory becomes exact (as linear theory gets 
more and more accurate), and 5(x,z) e(,g(x, 2) = D{z) eo(q), as 2 — > 00. Because of our normalization of D, here 
£o(q) is the mass density fiuctuation linearly extrapolated to the present time and filtered on the background scale 
Ro- The background smoothing scale Ro actually has a twofold role in our analysis. In the linear theory approach 
of Section 2 it was introduced and required to be much larger than the halo size, in order to get a self-consistent 
definition of halo counting fields, with the desirable feature of being free of the cloud-in-cloud problem. In the present 
non-linear analysis, however, the background mass scale must be chosen large enough to ensure that the halo velocity 
field coincides with the one of the matter. 

The Lagrangian density contrast of haloes identified by a PS-type algorithm can be obtained from eq.(|2^) as 
5h{q\M, Zf) = {q, z\M, Zf)/npg{z\M, Zf) — 1, which leads to 

-3/2 

Sh{ci\M,zf) = 



1 - 



exp 



eo(q) 



- 2eo(q)<5c/0(2/) + S'ia'jD{zf)\l 



For 2> (Jo this expression simplifies to 



5fe(q|M, Zf) 



Dizf)eo{ci) 



exp 



6o(q)2 -2eo(q)5c/D(2/) 



2a2 



E 



b^eiM,Zf) 



eo(q) 



The first four Lagrangian bias factors evaluated at 2 = read 



bi:,{M,Zf)=D{zf) 



<5c 



bUM,Zf) = 



Dizf)^al 



bUM,Zf) = 



bo4{M, Zf) 



D{zf) 



6 5c 



D{zf)^cjt, 
St 



D{zf)^al 
10 5^ 



+ 



4-15 



(42) 
(43) 

(44) 
(45) 
(46) 
(47) 



D{zf)^at, D{zf)^al 

Note that, in full generality, b^dM, Zf) = D(z)^ bf{z\M, Zf). Adding the further requirement that the local fluctuation 
on the background scale Ro has not collapsed yet by the time of halo formation would make our object number density 
semi-positive definite both at the Lagrangian and Eulerian level, i.e. 5h > —1, at any time, only provided eo < tf. 

The general expression for the Lagrangian halo density contrast of eq.(|42[) is plotted in Figure 2 as a function 
of the background density field, for different halo masses. In the high-mass case positive mass fluctuations typically 
correspond to positive values of the Lagrangian halo density contrast (and viceversa) , while the trend is the opposite 
at low masses. The transition, once again, corresponds to halo masses around M*, in which case positive values 
of 5h only occur in regions with background density close to the mean. Also shown are two approximations to the 
Lagrangian halo density contrast obtained by expanding eq.(^^ up to first and second order in the background field. 
Except for halo masses near Af, , where a quadratic bias is clearly needed, a linear Lagrangian bias generally provides 
an accurate fit to <5;i(q) within the bulk of the to distribution. 
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U, V. u„ 

Figure 2. The exact expression for the Lagrangian halo density contrast of eq.(42) (solid lines) is plotted as a function of 
i-'o = eo/uo. The three panels refer to values of the halo masses such that f^j/t^ = 1/4 (left), ""^j/*^ = 1 (centre) and = 4 

(right), with t/ = l.Gd/ D(zf). The background mass scale is chosen so that = 0.01 cr^ . Also plotted are two estimates of Sf^ 
obtained by expanding the RHS of eq.(42) up to first (dashed lines) and second order (dotted lines) in eo. Note that, because 
of our choice of variables, all the curves are independent both of Zf and f2o. 



The Eulerian bias field finally reads 
h'={^,z\M,Zf) = 1 + ^ t/^Y^ h^iq\M,Zf)eo(ci) . (48) 

It can be seen that, in the linear regime, where (S(x, z) « D{z) eo(q) ^ 1, the expression for in MW [i.e. our eq.(p6[)] 
is recovered, provided b^ (q|Af, z/) is replaced by its first-order approximation bai(AI, Zf). It is however important to 
realize that the exact expression in eq.(p8|) implies that the Eulerian bias field of dark matter haloes b^(x, A/, z/) 
is both non- linear, in that it depends on 5(x), and non-local, as it depends on the Lagrangian position q through 
bo (q|Af, Zf) eo(q), simply because of inertia. 

Our exact expression for the Eulerian halo bias, eq.(p8|), generally involves quantitative corrections to the MW 
approximate bias formula. In some cases, however, the MW relation may even fail to predict the correct qualitative 
behaviour of the halo-to-mass bias. This is the case, in fact, of those initially underdense fluid elements in Lagrangian 
space, £o(q) < 0, which, after an initial expansion phase, turn around to undergo a phase of local compression, 
so that the corresponding Eulerian fluid element eventually becomes overdense, (5(x(q, 2), z) > 0, and collapses. 
This is a well-known non-linear effect caused by the shear component of the velocity fleld, i.e. by the tidal force 
of the surrounding matter. For Gaussian initial conditions, the occurrence of such an event can be estimated by 
the Zel'dovich approximation as affecting 42% of the overall Lagrangian volume (Doroshkevich 1970; Shandarin & 
Zel'dovich 1984); Hui and Bertschinger (1996), using a different approximation, estimated this effect as affecting at 
least 39% of the total Lagrangian volume. In all such cases the MW formula would incorrectly predict bias instead 
of anti-bias, for halo masses M > Af», and anti-bias instead of bias, for Af < A/,. The problem may be generally 
less severe than the above heuristic argument would suggest, as, at a flxed epoch z, only a smaller fraction of such 
Lagrangian patches have already turned around from their initial expansion; this is even more true for large mass 
haloes, which probe the underlying mass distribution in a more linear regime, where the MW formula gets closer to 
the exact one. As a tentative conclusion let us say that one should be careful in applying the linear MW bias [i.e. our 
eq.(p^)] at the Eulerian level especially in connection with halo masses much smaller than Af«. 

The most important application of eq.(^8|) is that it allows to generate Eulerian maps of the local comoving halo 
number density per unit mass, npg(M, Zf)[l + 5h{^, z\M, zj)], given the non-linearly evolved mass density contrast 
5(x, z) (with Lagrangian resolution _Ro) and the corresponding Lagrangian mass and halo density fluctuation fields, 
eo(q) and 5h{q\M, Zf), respectively. 

In order to account for halo merging, at this level, one just has to assume a suitable link between the formation 
and observation epochs, which, in the simple PS theory amounts to the replacement Zf z, in the above expressions 
for b^ 
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Recalling that mass conservation can be recast in terms of the Jacobian determinant J = ISx/c'qll of the mapping 
q X, as 1 + <5[x(q, z), z] — J(q, z)~^ , one finds the exact relation 

b^ (x(q, z),z\M, Zf) = ! + [!- J(q, b^ (q|M, Zf) 6o(q) . (49) 

It can be useful to illustrate the meaning of this expression by considering various approximations to the evolution 
of the mass density in the non-linear regime, i.e. to the particle trajectories x(q, 2). Such approximation schemes 
should be thought of, not as self-consistent perturbative approaches to the actual dynamics, but as 'clever tricks' 
able to catch some aspects of the true dynamics, at least in the mildly non-linear regime. A detailed and systematic 
comparison of the performance of several approximations for different choices of the initial conditions has been made 
by Sathyaprakash et al. (1995). 



3.1.1 Zel'dovich approximation 

In the Zel'dovich approximation (ZEL; Zel'dovich 1970) the displacement vector is S = — D Vq(po(q), where 'fi{c\) is 
the linear peculiar gravitational potential, suitably rescaled so that Vqiy9o(q) = eo(q). Indicating by AQ,(q) (a — 1, 2, 3) 
the eigenvalues of the deformation tensor d^(po{ci)/dqa dq/s, we obtain for the Eulerian bias field 

.A N , eo(q)b^(q|M,2j) _^ , b^ (q| Af , 2/) n/..^M2(q) , n..^2 M3(q) ^ 



biBi(x(q, z), z\M, Zf)^l + ^7 r , .1 = 1 + 



.U{z)^^ + D{z 



/ii(q) /ii(q) 



(50) 



i-nLi[i--D(^)^"(q)] 

Here /ii(q) = Ai -I- A2 -I- A3 = eo(q), ^2 = A1A2 -I- A1A3 -I- A2A3 and ^3 = A1A2A3 are the three invariants of 
the deformation tensor. If one makes the further approximation of replacing the Lagrangian bias by its first-order 
estimate of eq.(^), it can be checked that the expression of hzEL coincides with the MW result, both at sufficiently 
early times {D ^ 1) and in the case of one- dimensional perturbations, for which fJ,2 ~ = fj,^ and the Zel'dovich 
approximation represents the exact solution to the non-linear dynamics. 

It is important to stress that we are not forced to take the above result as a perturbative expression. An accurate 
approximation to the Eulerian bias field would in fact consist in evolving the mass according to the truncated (on the 
scale Mo) Zel'dovich approximation (Kofman 1991; Kofman et al. 1992; Coles, Melott & Shandarin 1993) and using 
the full expression for the Lagrangian bias. Being a random field, the Eulerian halo bias is completely characterized by 
a probability density functional, thus for a given mass M and formation redshift Zf there exists a whole distribution of 
possible values of b^, related to the particular environment where the object forms as well as to the initial conditions 
leading to that site. Starting from the ZEL expression in eq.(|5o|) one could explicitly obtain the probability distribution 
function pih^Eh) dhzsL, by integrating over the joint distribution of the invariants Ha [an expression for the latter 
is given in Kofman et al. (1994)]. These specific applications of our results will be discussed elsewhere. 

Equation ( |5o| ) has the merit of clearly displaying the intrinsic non-locality of the Eulerian bias. Only in some 
simplified cases there exists a local mapping between b^ and 5, so that an expansion of the halo density contrast in a 
hierarchy of Eulerian bias factors, bf, bf, etc., makes sense. One example is provided by the linear theory approach 
of Section 2.6; further examples are given below. 



3.1.2 Frozen-flow approximation 

According to the frozen-flow approximation (FFA; Matarrese et al. 1992) the Eulerian density field can be written as 

l + <5(x(q,2),2) =exp / dD eo[:ic{ci, D)] , (51) 
Jo 

where the integral is calculated along the trajectory of the fluid element. Note that, since in the frozen-flow approx- 
imation shell-crossing never occurs, the mapping q ^ x can be inverted at any time. The solution ( |5l|) might be 
replaced in eq.(^8|) to obtain a non-local expression for the FEA bias parameter. However, we can make a further 
step by noting that, for Lagrangian points q^ corresponding to local extrema of the initial gravitational potential, 
Vq<^o(q, ) ~ 0, FFA predicts x, — x(q, ,2;) — q, , and 
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1 + <5(x,, z) = exp [l>(2;) eo(x,)] 



(52) 



One can speculate that such points represent the preferential sites for the formation of massive haloes, which could 
be associated to clusters of galaxies, and use this approximate expression to obtain 

l + 5(x.,2) , r . b^(x*|A/,zy) 



b^f^(x.,2|M, Zf)^l + - 



In [l + 5{y 



(53) 



<5(x.,^) L V " '\ D{z) 

Expanding this expression in powers of S, to first-order we recover the MW expression, eq.(|26|), while to second-order, 
we obtain 



{z\M,Zf) 



1 



D{zf) 
D(z) 



(54) 



which differs from the linear theory prediction of eq.(|2^). Analogous results could be obtained using the frozen- 
potential approximation (Brainerd, Scherrer & Villumsen, 1993; Bagla & Padmanabhan 1994), with the main differ- 
ence that the S evolution would be slowed down compared to FFA. Quite interesting is that the lognormal model by 
Coles and Jones (1991) assumes that the quantity 1 -f- (5(x, z) can be always approximated by the exponential of the 
linear density field at the same Eulerian position, so that the expressions above for the Eulerian bias factor in FFA 
would apply to all Eulerian points x. Of course, the validity of these approximate expressions for the bias should be 
checked against the results of N-body simulations. 



Another way to get a local mapping between the evolved halo density field and the underlying matter perturba- 
tions is to approximate the non-linear evolution of the mass by the spherical top-hat model. This method has been 
followed by MW and Mo et al. (1996). 



3.2 Perturbative evaluation of the Eulerian halo density contrast 

In the previous section we demonstrated that the Eulerian bias is a non-linear and non-local function of the density 
fluctuation fleld. The 'non- locality', in particular, comes from the fact that the halo number density fluctuation in x 
is determined by the initial halo number fluctuation at the Lagrangian position q, which, in turn, is related to the 
linear mass fluctuation in the same point, through a hierarchy of Lagrangian bias parameters. Here we want to derive 
an approximate expression for 5h(x(q, D), D), by applying the second-order Eulerian perturbation theory. Whenever 
it will be necessary to go from the Lagrangian position q to the Eulerian one, the Zel'dovich approximation will show 
sufficient. 

Within the linear regime, the Eulerian solution of the continuity equation is simply, 5^^\-k,D) — Z)eo(x). 
The mildly non-linear regime may be approximately described by the second-order solution (Bouchet et al. 1992; 
Bernardeau 1994; Catelan et al. 1995) 

5(^'(x,Z)) = i [l--^j5(i'(x,7?)^-7?uW(x).V5(i)(x,Z)) + iD^ + 9.4'^(x) 9„4^'(x) , (55) 



2 



D2 



2 

in such a way that S — 5'^' + (j'^' and higher-order corrections are neglected. Here u'^' (x) — —Vifo (x) is the (scaled) 
linear peculiar velocity, and ipo the (scaled) peculiar gravitational potential, linearly extrapolated to the present time. 
The second-order growth factor E — E[D) is quite a complicated function of D{Q.) (see Appendix A, for its explicit 
expression), but, in the vicinity of SI = 1 (actually in the range 0.05 < SI < 3), it can be approximated by the 
expression E ^ ^'in-'^/^^D'^ + 0[{n - 1)^] (see Bouchet et al. 1992). Therefore, the previous second-order solution 
is well approximated by the expression which holds in the Einstein de Sitter universe, namely (Fry 1984) 

6''^\yi,D) ^^5'-^'>{x,Df - Du'-^'>{x)-V5^''''{x,D) + ^D''d^u^^\x)d^uf\x) . (56) 

We want now to compute the corresponding second-order perturbative correction, 5[^''(x,_D), to the linear halo 
density fluctuation fleld, S'^\'k, D). From equation (^) we obtain 

S, ^ + Si'^ + fef 5« + + i 4''' , (57) 
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where, to maintain compact the notation we wrote e.g. = (5'-''(x, D). The Lagrangian bias factors bf — bi{z\M, Zf) 
and &2 — ''2(21 Af, z/) are those given in eqs.(|l^) and (|l7|). Notice that the perturbative expansion of 5^ holds at 
sufficiently early times and/or large scales, while the validity of the expansion of 5h(q) in powers of eo(q) is based on 
assuming a suitably large smoothing radius Ro on the background field e(q). 

The key point is that the first-order density field at the Lagrangian position q originates a non-local term, when 
written at the corresponding Eulerian position x. Using the Zel'dovich approximation x = q -f Du*^^', one obtains 



= Si^-* - Du^^^ ■ V5i^\ Finally, defining Su 



gets = (1 + hi) 5'^' and 



<5f = 



(1) 



2" \' • 752; """/S """^ • (58) 

Thus, the non-locality has the effect of modifying the inertia term u'^-* ■ V5'^', which gets multiplied by the factor 
(1 -I- 6^'). The dynamical properties of the random field 5h may be equivalently analyzed in terms of its Fourier 
transform Sh{]i.,t) where k is the comoving wavevector. Thus, the second-order solution may be written as 

cikidk2 



5f (k,Z))=75^ 



5„ (ki k2 - k) Hi,'' (ki , ka ; , 6^ ; n) <5i (ki ) 5i (k2 



(27r)3 



(59) 



where the symmetrized kernel Tig 



(2) 



reads 



H 



(2) 



(ki,k2;bf ,62 ,!^) 



hi 



fcl ^ fc2 

k2 ki 



The corresponding kernel for the Einstein-de Sitter case reads 



ki-ka 1 

fci k2 2 



1 + 



ki ks 



f- 

V ki 



k2 



(60) 



Hf\ki,-k2;btb^,n- 



1) 



^+b,+-b2 



+ 



1 + bi: / fci 



/fci ^ fca \ ki ■ ka ^ 2 / ki ■ ka \ 
V fc2 fci y fci fca 7 V fci fca / 



(61) 



3.3 Halo bispectrum and skewness 

A possible application of these results is the evaluation of the bispectrum and corresponding skewness of the halo 
distribution. A related calculation has been performed by Fry (1996), who assumed the bias to be local in Eulerian 
space at Zf. It should be clear that our model is quite different to the local Eulerian bias prescription applied to 
the analysis of the skewness by Fry and Gaztafiaga (1993). Moreover, the latter treatment, unlike ours, lacks of 
any prediction for the value of the difi'erent bias parameters. We recall that the value of the gravitationally induced 
skewness of the mass is 

5-^=4-2^, (62) 
for unfiltered fields, and 

^W = ^-4-2A_^(7^), (63) 

for a spherical top-hat filter, where 7 = —dlna{TZ)^ /dlnTZ (Bernardeau 1994). The smoothing radius TZ should not 
be confused with R, defining the halo mass: one is obviously interested in computing the skewness on a smoothing 
scale much larger than the typical size of the single objects. In the Einstein-de Sitter universe, and for a scale-free 
power-spectrum, with spectral index n, the latter reduces to S{TZ) — 34/7 — (n + 3), for —3 < n < 1. 

The derivation of the halo skewness {(Jj^) ~ •^{Sf^^^Sl^^) is simple. Assuming that the Eulerian halo density field 
is smoothed by a top-hat filter, the halo skewness parameter Sh is, for a generic value of Q, 

^ t-v _ _ 4-2^+6 b'^{z\M, zj) + 3 bi{z\M, Zf) - [1 + bi{z\M, zj)] 7(7^) 

' ^ 1^ [i+biiz\M,z,)r • ^''^ 

The asymptotic value of Shiji-', z, Q), for a fixed formation redshift Zf, is S — 7(7?.) as z ^ —1, in the open and flat 
cases, while, for f2 > 1, this value is attained at the time of maximum expansion, corresponding to z — —1/ilo. This 
limit gives the value of the underlying mass skewness: in the absence of merging the haloes would evolve towards an 
unbiased distribution in the far future. It is of interest to write the halo skewness in the Einstein-de Sitter universe 
and for a scale-free linear power-spectrum. 
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f+6bnz\M,Zf) + 3bUz\M,Zf)-{n + 3)[l + bHz\M,Zf)] 
Sh{n;z,n=l) = 1 L_ (65) 

[l + b^{z\Al,Zf)] 

As for the mass skewness, the dependence on the smoothing scale TZ now simply translates into a dependence on the 
spectral index n. Again, the standard value for the mass skewness 34/7 — (n-|- 3) is always recovered at the end of the 
expansion phase. The skewness parameter is shown in Figure 3 for different values of fJo and for a scale-free model 
with n = —2. For objects observed at the present time, « = 0, wo vary the collapse epoch Zf, which may simulate 
different models of gala^xy formation inside dark haloes. By varying together z = Zf, we instead show the skewness 
evolution in the instantaneous merging model. We also consider the case of varying only z: this gives the evolution of 
the skewness in a model in which the objects did not suffer amy merging after their formation at Zf. Finally, we show 
the evolution of the skewness parameter of filtered mass fluctuations; note that the Einstein-de Sitter case displays no 
redshift dependence, simply because of self-similarity; for sensible values of fio 7^ 1 also the mass skewness of non-flat 
Friedmann models experiences very little evolution. The redshift dependence of Su is therefore mostly due to that of 
the Lagrangian bias factors. Quite interesting, in this respect, is the fact that the halo skewness plotted in the two 
top panels of Figure 3 displays a turning point in its redshift dependence: this typically occurs when M w Mt{zf). 
Of particular interest is also the expression for the halo bispectrum defined by the relation 

(5h{k2, D) 5h(k2, D) 4(k3, D)) = {27^)^5^ (ki + ka + kg) Bh(ki, ka, kg; . (66) 

The leading term shows the characteristic hierarchical pattern 

Bh(ki,k2,k3;£)) [l + hi{z\M,Zf)Y [2 7Y^'^(ki, ka; fef , 6^, 0) P(A;i) P(fc2) + cyclic terms] , (67) 

where P{k) is the primordial density power-spectrum defined by {5i (fci ) 5i (^2 )) — {2-Ky''5^ (ki -|- k2) P(/ci ) , and the 
two cyclic terms are obtained by the substitutions ki — » k2 , ki — » ka and k2 — » ka . Typically, as for the hierarchical 
mass bispectrum, the halo bispectrum is largely scale dependent, while its dependence on the ft-shape is rather weak. 
One way to eliminate the scale dependence and look at the residual shape dependence is to analyze the 'effective' 
bispectrum amplitude Q (Fry 1984) 

^ Ph{ki,D)Ph{k2,D) + Ph{ki,D)Ph{k2,D) + Ph{k2,D)Ph{k3,D) - ^ ' 

The halo power-spectrum is biased with respect to the mass one, Ph{k,D) = D'^ [l -|- &{'(z)] ^ P(fc). For a power- 
law spectrum, the amplitude Q generally depends on the spectral index n, owing to the wavenumber modulation 
introduced by the kernel H^g\ki,k.2) (cf. Figure 4). For equilateral triangle configurations, Q gets an n-independent 
value, namely 

O (0. - Hi - 3 + 2 b'^iAM, zj) + b^{z\M, Zf) 
reducing to 

Qeq{0.= l-Z) = -2 , (70) 

[l + bi{z\M,Zf)] 
in the Einstein-de Sitter universe. 



3.4 Local Lagrangian bias 

So far, our model has been treated as being fully predictive. Once the cosmological scenario and the structure formation 
model have been fixed, our algorithm contains no fitting parameters. This is because we used a local version of the 
PS theory to generate the Lagrangian halo density contrast. One could, however, take a more general point of view 
and assume that the Lagrangian halo density contrast 5h (q) is specified in terms of the linear background density 
field e6g(q, z) = L'(z)eo(q) by a set of unknown bias parameters bf{z), as follows, 

e=i ' e=i 
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Figure 3. The filtered skewness parameter is plotted, for Qo = 0.5, 1, 1.5, for a scale-free model with n = —2. The halo masses 
are selected with the same linear mass variance = 10, corresponding to the same present-day bias parameters. We take 
everywhere 5c = 1.69. The top left panel refers to objects observed at z = 0, with varying formation redshift Zf. The top right 
panel shows the effect of varying simultaneously z = Zf. In the bottom left panel we fix Zf = 5 and look at different observation 
redshifts z < Zf. The bottom right panel, finally, shows the evolution of the skewness parameter of filtered mass fluctuations. 



Defining now bi = bi{z) = 1 + (z) and 62 = ^2(2) ~ 2bi{z) + 62 (z), according to eq.(|2q), and replacing these 
expansions in our previous treatment, we recover the general expression ( |59| ) for the second-order halo density contrast, 
with the more general kernel 



H\:'{ki,k2;bi,b2,n) = - 



2 U2 ki 



ki ■ k2 



1 E 



ki ■ k2 



(72) 



fci fc2 2 V D V V fci fc2 

Comparing this relation with the analogous one obtained with a local Eulerian bias expansion (e.g. Fry, Melott & 
Shandarin 1995; Matarrese, Verde & Heavens 1997), we see that the bispectrum for a set of objects selected by a 
local Lagrangian bias differs from the results of the local Eulerian bias by the extra inertia term 
61 — 1 / fci , ^2 \ ki • k2 



2 U2 fci. 



fci k2 



(73) 



which implies a different shape dependence. 



© 0000 RAS, MNRAS 000, 000-000 



20 P. Catelan, F. Lucchin, S. Matarrese and C. Porciani 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 



Figure 4. The halo bispectrum amphtude Q{9) for configurations with sides ki = I, k2 = 1/2, separated by an angle 6 is 
plotted vs. 6 for scale-free models with n = — 2 and n = —1 at z = z f = and in a fiat Universe. Two cases are shown for each 
panel: the local Lagrangian bias model, with linear Eulerian parameters bi = 2 and 62 = 1 {solid line) and the local Eulerian 
bias model, with the same bias parameters {dashed line). 



The halo bispectrum amphtude Q{6), nt z = Zf = 0, for configurations with sides ki = 1, k2 = 1/2, separated 
by an angle 6, is shown in Figure 4, for scale- free models with n — ~2 and n = — 1, with f2 = 1. Two different cases 
are considered: our local Lagrangian bias model, with linear Eulerian parameters bi — 2 and 62 = 1 and the local 
Eulerian bias model of Fry and Gaztafiaga, with the same Eulerian bias parameters. 

Similar reasoning would apply to the skewness, for which the local Lagrangian vs. Eulerian bias hypothesis 
implies a change of the scale dependence, through the extra term 

With adequate modeling of galaxy formation inside dark matter haloes (e.g. Kauffmann et al. 1997, and references 
therein) the results of this section can be used to predict the clustering properties of galaxies at different redshifts. In 
particular, the specific shape dependence of the bispectrum (and related scale dependence of the skewness), implied 
by our local Lagrangian bias prescription, would reflect into a detectable signature in the statistical properties of 
the galaxy distribution. Our model, therefore, provides a valid alternative to local Eulerian bias schemes (e.g. Cen & 
Ostriker 1992; Coles 1993; Fry & Gaztafiaga 1993; Catelan et al. 1994; Mann et al. 1997). 



4 CONCLUSIONS 

In this paper we studied the non-linear evolution of the clustering of dark matter haloes, using a stochastic approach 
to single out the halo formation sites directly in Lagrangian space. Our model is based on a local version of the 
Press- Schechter theory, which becomes free of the cloud-in-cloud problem after a suitable coarse-graining procedure 
is applied. The non-linear evolution of the halo distribution is then followed exactly by relating it to the dynamics of 
the Lagrangian patch of fluid which the nascent halo belongs to. 

This formalism allowed us to obtain the bias random field relating the local halo density contrast to the underlying 
mass distribution. The expression for the halo bias field, reported in eqs.(^) and (^9|), represents the most relevant 
result of our paper. Because of the locality in Lagrangian space inherent in our approach, such a bias field turns 



© 0000 RAS, MNRAS 000, 000-000 



Halo bias field 21 



out to be non-local in Eulerian coordinates, which has relevant implications for the clustering properties of luminous 
objects like galaxies and galaxy clusters that formed inside dark matter haloes. 

Our method contains two Lagrangian smoothing scales. The scale R, selecting the halo mass, and the background 
scale -Ro 3> -R allowing us to define the Lagrangian halo counting field as the local PS mass function in a patch with 
comoving background density Qb [1 + £69 (q, z)], ttg being the linear mass fluctuation smoothed on the scale Ro. Given 
the role of the latter, it would appear that our description of halo clustering makes sense only on scales larger than 
Ra- On the other hand, the derivation of the Lagrangian correlation function in Section 2.4, which does not make 
use of the background field, suggests that we can actually extrapolate our Lagrangian results down to separation 
comparable to the halo size. This result is further confirmed by an analysis in terms of space-correlated Langevin 
equations (Porciani et al. 1997). The numerical results of MW and Mo et al. (1996) support the idea that such 
an extrapolation would apply even in the non-linearly evolved case. In our treatment of the non-linear regime, the 
background scale Ro appears with a complementary role. It is the minimum scale ensuring that the nascent haloes 
are indeed comoving with the Lagrangian fluid patch which they belong to. This would reasonably require that the 
Lagrangian fluid elements evolve with negligible orbit crossing (e.g. Kofman et al. 1994). 

Once again, let us stress that our approach makes no assumptions about the merger rates of the considered 
objects. The clear distinction between observation and formation redshift, z and z/, in our approach implies that the 
instantaneous merging hypothesis, implicit in the standard PS model, as well as any other realistic approximation 
can be easily accommodated into our scheme as just the way to relate zj and z. 

Our method for evolving the spatial distribution of the haloes is indeed much more general than the specific 
application we have considered so far. Given any Lagrangian population of objects specified by some set of physical 
properties M (like mass and formation threshold in our halo model) , with conserved mean comoving number density 
fiobjiM) and local Lagrangian density contrast 5obj(q|A4), our results imply that, at any redshift z, their comoving 
local density in Eulerian space is given by 



where x(q, z) = (\ + S(q, z), and S(q, z) is the displacement vector of the q-th Lagrangian element. Smoothing the 
initial gravitational potential on some scale -Ro is again required, so that the objects assigned to the q-th patch 
can be sensibly assumed to be comoving with it. This method could be used, for instance, to follow the clustering 
of the Lagrangian density maxima in the non-linearly evolved mass density field. This suggests that, starting from 
low-resolution numerical simulations, one can generate mock catalogues of the given class of objects, with local 
density correctly specified up to some resolution scale. One can understand the last relation as a local version of the 
Chapman-Kolmogorov equation of stochastic processes (e.g. van Kampen 1992), stating that the local Eulerian object 
distribution is the convolution of the Lagrangian object density with the 'conditional particle density', (5_d[x — x(q, z)\, 
i.e. the probability that a particle is found in x at redshift z given that it was in q as z — > 00. The only underlying 
hypothesis being, once again, that these objects move exclusively by the action of gravity. It may be worth to notice 
that the latter equation is actually more general than eq. (|38|), as it also holds in the presence of multi-streaming. 

Our non- linear stochastic approach can be already considered successful in that, besides recovering the PS mass 
function, it provides a self-consistent derivation of the Eulerian halo bias, which, to a first approximation, reduces 
to the MW formula. We, however, also predict both quantitative and qualititive corrections to the MW results, that 
clearly need to be checked against the outputs of numerical simulations. A definite prediction of our analysis is, for 
instance, the form of the skewness and of the bispectrum of the spatial halo distribution, which significantly deviates 
from that deduced with any local Eulerian bias prescription. 
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APPENDIX A: GROWTH FACTORS IN FRIEDMANN UNIVERSE MODELS 



The expressions for the first and second-order growth factors D{z) and E(z) have not been given in the main text. 
An easy derivation can be given following Shandarin (1980) and using the relation 



- 1 = {n-' - i) {1 + z)-' . 



(Al) 



We consider only cases with vanishing cosmological constant. The growth factor 'D(«; Ho) of linear density perturbar- 
tions reads, for the different geometries, 



5 , 15 Slo(l+z) 

2 ' 2 1-flo 

_5 I 15 no(l-Hz) 
2 2 Qo-l 



Jloil + Z) 



t)-2,y(l-Uo)(i+noz) , 



1 + 



arctan — 



l + fioz 



(fio < 1) 
(fio = 1) 
(fio > 1) . 



(A2) 



The expressions for the second-order growth factors £{z;flo) are slightly more cumbersome 



£{z;Qo)= - 



25 225 no{l + z) 
y ~ ~8 l-Oo 
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' no(i-f2) 
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(Slo < 1) , (A3) 
(A4) 
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(no > 1) 



(A5) 



Notice that we are implicitly adopting here the normalization suggested by Shandarin (1980), so that, in the 
limit — > oo one recovers the Einstein de Sitter case, 'D(2;no) — » (1 -|- z)~^. However, in the main text we 
normalized to unity the linear growing factors extrapolated to the present time; so, for any geometry, we define 
D{z) =V{z;no)/'D{z = 0;flo) and E{z) = £{z■,flo)/[V{z = 0■,no)f■ 
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